home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_y.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  7.7 KB  |  289 lines
  1. /* specfunc/bessel_y.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_pow_int.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_trig.h>
  28. #include <gsl/gsl_sf_bessel.h>
  29.  
  30. #include "error.h"
  31.  
  32. #include "bessel.h"
  33. #include "bessel_olver.h"
  34.  
  35. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  36.  
  37. /* [Abramowitz+Stegun, 10.1.3]
  38.  * with lmax=15, precision ~ 15D for x < 3
  39.  *
  40.  * checked OK [GJ] Wed May 13 15:41:25 MDT 1998 
  41.  */
  42. static int bessel_yl_small_x(int l, const double x, gsl_sf_result * result)
  43. {
  44.   gsl_sf_result num_fact;
  45.   double den = gsl_sf_pow_int(x, l+1);
  46.   int stat_df = gsl_sf_doublefact_e(2*l-1, &num_fact);
  47.  
  48.   if(stat_df != GSL_SUCCESS || den == 0.0) {
  49.     OVERFLOW_ERROR(result);
  50.   }
  51.   else {
  52.     const int lmax = 200;
  53.     double t = -0.5*x*x;
  54.     double sum = 1.0;
  55.     double t_coeff = 1.0;
  56.     double t_power = 1.0;
  57.     double delta;
  58.     int i;
  59.     for(i=1; i<=lmax; i++) {
  60.       t_coeff /= i*(2*(i-l) - 1);
  61.       t_power *= t;
  62.       delta = t_power*t_coeff;
  63.       sum += delta;
  64.       if(fabs(delta/sum) < 0.5*GSL_DBL_EPSILON) break;
  65.     }
  66.     result->val = -num_fact.val/den * sum;
  67.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  68.  
  69.     return GSL_SUCCESS;
  70.   }
  71. }
  72.  
  73.  
  74. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  75.  
  76.  
  77. int gsl_sf_bessel_y0_e(const double x, gsl_sf_result * result)
  78. {
  79.   /* CHECK_POINTER(result) */
  80.  
  81.   if(x <= 0.0) {
  82.     DOMAIN_ERROR(result);
  83.   }
  84.   else if(1.0/GSL_DBL_MAX > 0.0 && x < 1.0/GSL_DBL_MAX) {
  85.     OVERFLOW_ERROR(result);
  86.   }
  87.   else {
  88.     gsl_sf_result cos_result;
  89.     const int stat = gsl_sf_cos_e(x, &cos_result);
  90.     result->val  = -cos_result.val/x;
  91.     result->err  = fabs(cos_result.err/x);
  92.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  93.     return stat;
  94.   }
  95. }
  96.  
  97.  
  98. int gsl_sf_bessel_y1_e(const double x, gsl_sf_result * result)
  99. {
  100.   /* CHECK_POINTER(result) */
  101.  
  102.   if(x <= 0.0) {
  103.     DOMAIN_ERROR(result);
  104.   }
  105.   else if(x < 1.0/GSL_SQRT_DBL_MAX) {
  106.     OVERFLOW_ERROR(result);
  107.   }
  108.   else if(x < 0.25) {
  109.     const double y = x*x;
  110.     const double c1 =  1.0/2.0;
  111.     const double c2 = -1.0/8.0;
  112.     const double c3 =  1.0/144.0;
  113.     const double c4 = -1.0/5760.0;
  114.     const double c5 =  1.0/403200.0;
  115.     const double c6 = -1.0/43545600.0;
  116.     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*c6)))));
  117.     result->val = -sum/y;
  118.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  119.     return GSL_SUCCESS;
  120.   }
  121.   else {
  122.     gsl_sf_result cos_result;
  123.     gsl_sf_result sin_result;
  124.     const int stat_cos = gsl_sf_cos_e(x, &cos_result);
  125.     const int stat_sin = gsl_sf_sin_e(x, &sin_result);
  126.     const double cx = cos_result.val;
  127.     const double sx = sin_result.val;
  128.     result->val  = -(cx/x + sx)/x;
  129.     result->err  = (fabs(cos_result.err/x) + sin_result.err)/fabs(x);
  130.     result->err += GSL_DBL_EPSILON * (fabs(sx/x) + fabs(cx/(x*x)));
  131.     return GSL_ERROR_SELECT_2(stat_cos, stat_sin);
  132.   }
  133. }
  134.  
  135.  
  136. int gsl_sf_bessel_y2_e(const double x, gsl_sf_result * result)
  137. {
  138.   /* CHECK_POINTER(result) */
  139.  
  140.   if(x <= 0.0) {
  141.     DOMAIN_ERROR(result);
  142.   }
  143.   else if(x < 1.0/GSL_ROOT3_DBL_MAX) {
  144.     OVERFLOW_ERROR(result);
  145.   }
  146.   else if(x < 0.5) {
  147.     const double y = x*x;
  148.     const double c1 =  1.0/6.0;
  149.     const double c2 =  1.0/24.0;
  150.     const double c3 = -1.0/144.0;
  151.     const double c4 =  1.0/3456.0;
  152.     const double c5 = -1.0/172800.0;
  153.     const double c6 =  1.0/14515200.0;
  154.     const double c7 = -1.0/1828915200.0;
  155.     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7))))));
  156.     result->val = -3.0/(x*x*x) * sum;
  157.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  158.     return GSL_SUCCESS;
  159.   }
  160.   else {
  161.     gsl_sf_result cos_result;
  162.     gsl_sf_result sin_result;
  163.     const int stat_cos = gsl_sf_cos_e(x, &cos_result);
  164.     const int stat_sin = gsl_sf_sin_e(x, &sin_result);
  165.     const double sx = sin_result.val;
  166.     const double cx = cos_result.val;
  167.     const double a  = 3.0/(x*x);
  168.     result->val  = (1.0 - a)/x * cx - a * sx;
  169.     result->err  = cos_result.err * fabs((1.0 - a)/x) + sin_result.err * fabs(a);
  170.     result->err += GSL_DBL_EPSILON * (fabs(cx/x) + fabs(sx/(x*x)));
  171.     return GSL_ERROR_SELECT_2(stat_cos, stat_sin);
  172.   }
  173. }
  174.  
  175.  
  176. int gsl_sf_bessel_yl_e(int l, const double x, gsl_sf_result * result)
  177. {
  178.   /* CHECK_POINTER(result) */
  179.  
  180.   if(l < 0 || x <= 0.0) {
  181.     DOMAIN_ERROR(result);
  182.   }
  183.   else if(l == 0) {
  184.     return gsl_sf_bessel_y0_e(x, result);
  185.   }
  186.   else if(l == 1) {
  187.     return gsl_sf_bessel_y1_e(x, result);
  188.   }
  189.   else if(l == 2) {
  190.     return gsl_sf_bessel_y2_e(x, result);
  191.   }
  192.   else if(x < 3.0) {
  193.     return bessel_yl_small_x(l, x, result);
  194.   }
  195.   else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1.0)) {
  196.     int status = gsl_sf_bessel_Ynu_asympx_e(l + 0.5, x, result);
  197.     double pre = sqrt((0.5*M_PI)/x);
  198.     result->val *= pre;
  199.     result->err *= pre;
  200.     return status;
  201.   }
  202.   else if(l > 40) {
  203.     int status = gsl_sf_bessel_Ynu_asymp_Olver_e(l + 0.5, x, result);
  204.     double pre = sqrt((0.5*M_PI)/x);
  205.     result->val *= pre;
  206.     result->err *= pre;
  207.     return status;
  208.   }
  209.   else {
  210.     /* recurse upward */
  211.     gsl_sf_result r_by;
  212.     gsl_sf_result r_bym;
  213.     int stat_1 = gsl_sf_bessel_y1_e(x, &r_by);
  214.     int stat_0 = gsl_sf_bessel_y0_e(x, &r_bym);
  215.     double bym = r_bym.val;
  216.     double by  = r_by.val;
  217.     double byp;
  218.     int j;
  219.     for(j=1; j<l; j++) { 
  220.       byp = (2*j+1)/x*by - bym;
  221.       bym = by;
  222.       by  = byp;
  223.     }
  224.     result->val = by;
  225.     result->err = fabs(result->val) * (GSL_DBL_EPSILON + fabs(r_by.err/r_by.val) + fabs(r_bym.err/r_bym.val));
  226.  
  227.     return GSL_ERROR_SELECT_2(stat_1, stat_0);
  228.   }
  229. }
  230.  
  231.  
  232. int gsl_sf_bessel_yl_array(const int lmax, const double x, double * result_array)
  233. {
  234.   /* CHECK_POINTER(result_array) */
  235.  
  236.   if(lmax < 1 || x <= 0.0) {
  237.     GSL_ERROR ("error", GSL_EDOM);
  238.   }
  239.   else {
  240.     gsl_sf_result r_yell;
  241.     gsl_sf_result r_yellm1;
  242.     int stat_1 = gsl_sf_bessel_y1_e(x, &r_yell);
  243.     int stat_0 = gsl_sf_bessel_y0_e(x, &r_yellm1);
  244.     double yellp1;
  245.     double yell   = r_yell.val;
  246.     double yellm1 = r_yellm1.val;
  247.     int ell;
  248.  
  249.     result_array[0] = yellm1;
  250.     result_array[1] = yell;
  251.  
  252.     for(ell = 1; ell < lmax; ell++) {
  253.       yellp1 = (2*ell+1)/x * yell - yellm1;
  254.       result_array[ell+1] = yellp1;
  255.       yellm1 = yell;
  256.       yell   = yellp1;
  257.     }
  258.  
  259.     return GSL_ERROR_SELECT_2(stat_0, stat_1);
  260.   }
  261. }
  262.  
  263.  
  264. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  265.  
  266. #include "eval.h"
  267.  
  268. double gsl_sf_bessel_y0(const double x)
  269. {
  270.   EVAL_RESULT(gsl_sf_bessel_y0_e(x, &result));
  271. }
  272.  
  273. double gsl_sf_bessel_y1(const double x)
  274. {
  275.   EVAL_RESULT(gsl_sf_bessel_y1_e(x, &result));
  276. }
  277.  
  278. double gsl_sf_bessel_y2(const double x)
  279. {
  280.   EVAL_RESULT(gsl_sf_bessel_y2_e(x, &result));
  281. }
  282.  
  283. double gsl_sf_bessel_yl(const int l, const double x)
  284. {
  285.   EVAL_RESULT(gsl_sf_bessel_yl_e(l, x, &result));
  286. }
  287.  
  288.  
  289.